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ABSTRACT 



o 



We present a new numerical code that solves the general relativistic magneto-hydrodynamical (GRMHD) equations 
coupled to the Einstein equations for the evolution of a dynamical spacetime within a conformally-flat approximation. 
. This code has been developed with the main objective of studying astrophysical scenarios in which both, high magnetic 

' fields and strong gravitational fields appear, such as the magneto-rotational collapse of stellar cores, the collapsar model 

of GRBs, and the evolution of neutron stars. The code is based on an existing and thoroughly tested purely hydrody- 
namical code and on its extension to accommodate weakly magnetized fluids (passive magnetic-field approximation). 
■ These codes have been applied in the past to simulate the aforementioned scenarios with increasing levels of sophisti- 

Oh' cation in the input physics. The numerical code we present here is based on high-resolution shock-capturing schemes 

to solve the GRMHD equations, which are cast in first-order, flux-conservative hyperbolic form, together with the flux 
constraint transport method to ensure the solenoidal condition of the magnetic field. Since the astrophysical applica- 
tions envisaged do not deviate significantly from spherical symmetry, the conformal flatness condition approximation 
52 ' is used for the formulation of the Einstein equations; this has repeatedly shown to yield very good agreement with 

full general relativistic simulations of core-collapse supernovae and the evolution of isolated neutron stars. In addition, 
the code can handle several equations of state, from simple analytical expressions to microphysical tabulated ones. In 
^SJ ■ this paper we present stringent tests of our new GRMHD numerical code, which show its ability to handle all aspects 

^ ' appearing in the astrophysical scenarios for which the code is intended, namely relativistic shocks, highly magnetized 

' fluids, and equilibrium conflgurations of magnetized neutron stars. As an application, magneto-rotational core-collapse 

, simulations of a realistic progenitor are presented and the results compared with our previous flndings in the passive 

1/^ ■ magnetic-field approximation. 

o 

52 1. Introduction (see lKellv et al.|[2007i and references therein). In the case of 
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the mergers of two neutron stars, beheved to be the stan- 
The cohapse of rotatmg stellar cores and the mergmg of ^^^^ mech anism to account for hard-short GRBs, it was 
compact binaries (either neutron star-neutron star or neu- suggested (iPrice fc Rosswod IMoD that the strong shear 
tron star-black hole binaries) are two of the most important between the two neutron stars could generate strong mag- 



O 

> 

X, 

H , astrophysical scenarios involving compact objects, whose jj^Hq fields too 
. - - ' modeling requires the study of the dynamical evolution 

of a magnetized fluid in general relativity. In the case of A considerable effort has been made to develop spe- 

the rotational collapse of massive stellar cores, the mag- cial rela tivistic magneto - hydro d ynamics (SRM HD) codes 

netic field is thought to grow through the extraction of en- (see e.g.|Marti fc Miilleij (|2002D, |lbafiez| (|2006|) and refer- 

ergy from the diffe rential rotation generated during collapse ences therein) . Most works have considered the case of ideal 

(jMeier et al.lll976D . This idea is supported by the observa- MHD where the fluid is assumed to be a perfect conductor, 

tional fact that some neutron stars (magnetars) possess ex- In this case, the resulting system of equations is simplified 

tremely large magnetic fields (lO^^'-lOi^ Gauss), as inferred significantly, and can be solved by numerical codes designed 

from studies o f anomalous X-ray pulsa rs and soft gamma- specifically for hyperbolic syste ms. These codes iii clude the 

ray repeaters (jKouveliotou et al.l(l998l ). Furthermore, the use of Godunov-type schemes (|Komissaroy||l999|), numer- 

class of soft-long gamma-ray bursts (GRB) are probably ical techniques to keep the magnetic field divergence-free 

the result of jets formed in a rotational core-collapse event (see |T6th| [2000, and references therein), and efficient re- 

leading to a black ho le, according to the collapsar scenario covery schemes to derive primitive quantities from the con- 

(jWoosleY et al.lll993f) . In this case, the magnetic field most served ones (see [Noble etaD|2006|, and references therein), 

likely plays a crucial role in the formation and coUimation There is also a major activity in the development of codes 

of the jet. This scenario is supported by the existing correla- capable of simulating magnetized astrophysical flows in 

tion of some long GRB with core-collapse supernova events general relativity These codes integrate the ideal GRMHD 

equations for fixed background spacetimes using high-order 

Send offprint requests to: Pablo Cerda-Duran, conservative schemes based on either approximate or full 

e-mail: cerda@mpa-garching.mpg.de wave-decomposition Riemann solvers (jGammie et al.ll2003l : 
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Komissarovl 12005: Anninos et all l2005t lAnton et al.l [200l 



Del Zanna et al.ii2007tlTchekhovsko v et al."2007^. Resistive 
MHD flows were considered by [Komissarov (2007), and 
nonconservative GRMHD schemes and schemes relying 
on ar tificial viscosity were used by iDe Villiers fc HawlevI 
(|2003D and lAnninos et"al] (|2005D . Most codes have been 
applied to study disk accretion onto black holes and jet for- 
mation, but since the self-gravity of the fluid was not taken 
into account, these codes cannot simulate consistently the 
formation of the black hole and the evolution of the sur- 
rounding disk or torus. 

Only very recently, GRMHD codes are able to 
follow t he evolution of a dynamical spacetirn e . Th e 
code s of iDuez et al.l (l2005l). IShibata k SekiguchU (|2005D . 
and iGiacomazzo fc Rezzollal ( 20071) are based on the 
BSSN for mulation of the Einste i n equations for the 
spacetime (iNakamura et al.l 119871 : IShibata fc Nakamural 
119951 : iBaumgarte &: Shapiralll999I K high resolution shock- 
capturing schemes for the GRMHD equations involving ap- 
proximate Riemann solvers such as HLL or high-order cen- 
tral schemes, and on the constrai nt transport sch e me fo r 
the magnetic field. In the code of lAnderson et all (|2008f) . 
the Einstein equations are cast in first-order symmetric hy- 
perbolic form, and are solved using the generalized har- 
monic decomposition. While the code relies on the same 
type of GR hydrodynamics solvers as previously developed 
codes, it guarantees a divergence-free magnetic field by 
means of either projection methods or hyperbolic diver- 
gence cleaning. Both methods are easier to implement for 
the non-structured AMR grids employed in this code than 
the constrained transport method. All four codes rely on 
Cartesian co ordinates for three dim ensional simulations. 
The c odes of iDuez et al.l (|2005D and IShibata fc Sekiguchil 
(|2005f ) also provide the possibility to impose axisymmetry 
by means of the cartoon method for the spacetime evolu- 
tion and the use of cylindrical coordinates for the GRMHD 
equations. The equations of state (EOS) implemented in 
these codes consist of simple analytic expressions: poly- 
tropic EOS, ideal gas or hybrid EOS (see Sect. 12. 3p . One of 
the co des was extended to handle a tabulated microphysical 
EOS (IShibata et al.ll2007D . 

We present a new axisymmetric numerical code, cable of 
handling ideal MHD flows in dynamical spacetimes in gen- 
eral relativity, and designed particularly to investigate grav- 
itational core collapse. We use similar numerical schemes 
as in most of the other existing GRMHD codes (HRSC 
schemes and constraint transport), but we follow a simpler 
approach for the spacetime evolution. 

The new c ode is based on the hydrod ynamics code 
described in iDimmelmeier et"al] (l2002allbf) . and on 



its extensions discusse d in 
Cerda-Duran & Font. (|2007D . 
(|2007l) . The Maxwell equations 



Cerda-Duran et al.l (l2005t) 
and ICerda-Duran et all 



are ahead 



iaa y mo 
3 (120071) 



mcorpo- 
and 



rated in the code s of ICerda-Duran fc Font! 
ICerda-Duran et all (|2007f ). but only in the passive 
magnetic-field approximation, i.e. the contribution of 
the magnetic field to the energy-momentum tensor is 
neglected, and therefore has no impact on the dynamics. 
In the new code, we relax this assumption and incorporate 
magnetic field effects on the spacetime dynamics and 
the sel f-gravity of t h e flui d following the approach laid 
out in lAnton et al.l (|2006( ). The Einstein equations are 
formulated using the conformal flatness condition (CFG 
hereafter). This approximate treatment of the metric 



equations was first introduced by llsenberd (|1978D and 
I Wilson e t al.' (199(f), and was used to study rotational core 
collaps e (D immclmeier et al. 2002a,b), and binary neutron 
(lOec 



stars (|Oechslin etraLlT2007t ). Simulations with a second 
post-Newtonian extension of the CFG metric (named 
CFC-I-) showed small quantitative differences in the dy- 
namics and the gravitational waveforms (< 1%) compared 
to the CFC metric, both for rotational core collapse and 
for simulations of the e volution of single neutron stars 
(|Cerda-Duran et al.|[2005l) . Direct comparisons of the CFC 
approach wit h full general relati vi stic s imulat ions were 
reported by 'Shibat a fc Sekiguchil (|2004D and lOtt et al.l 
(2007a . b), who found that the differences in the collapse 
dynamics and the waveforms are minute demonstrating the 
suitability of CFC for performing accurate core collapse 
simulations. 

The CFC approach has some advantages compared with 
the BSSN formulation: (i) the Hamiltonian and momentum 
constraints of the Einstein equations are automatically sat- 
isfied, and (ii) the time step is less restrictive since it is 
determined by the largest fiuid eigenvalue, while in hyper- 
bolic formulations (such as BSSN) the largest eigenvalue is 
the speed of light. Consequently, the time steps are typ- 
ically ten times larger than those admitted in BSSN for 
the same grid for the evolution of neutron stars, and even 
larger during core collapse. However, there are also some 
disadvantages of the CFC approach. It neglects the gravi- 
tational wave content of the spacetime, i.e. when the grav- 
itational wave back reaction is important (e.g. in neutron- 
star mergers) the dynamics can not be modeled accur ately 
(e.g. compare the sim ulations of lOechslin et al.l (|2007( ) and 
IShibata fc Taniguchil (2006)). Furthermore, to compute the 
gravitational waveforms one needs to resort to the Einstein 
quadrupole formula, and, since the CFC metric equations 
are elliptic, the parallelization of the code for a large num- 
ber of processors is more difficult than in hyperbolic formu- 
lations such as BSSN, but still possible. 

Our code uses spherical polar coordinates and 
(presently) assumes axisymmetry. The most important ad- 
vantage of these coordinates with respect to Cartesian or 
cylindrical coordinates adopted by other numerical codes, 
is that they are more readily adapted to the astrophysical 
scenarios that we wish to study. Furthermore, it allows us to 
easily and properly cover the length-scales of core collapse 
ranging from the radius of the initial iron core (~ 1000 km) 
down to the radius of the neutron star (^ 10 km) by means 
of a logarithmically spaced radial grid. A disadvantage of 
our coordinate system concerns its possible extension to 3D 
because of the coordinate singularities at the center and at 
the axis. Moreover, as the azimuthal grid spacing decreases 
quadratically towards the axis with increasing grid reso- 
lution (for an equally-spaced angular grid), the Courant 
condition for the time step can be rather restrictive in 3D 
simula tions (a po s sible solution to this issue can be found 
e.g. in lZink et all (|2008f) ). 

The code can handle various equations of state rang- 
ing from simple analytical expressions (polytropes, ideal 
gas and hybrid EOS) to tabulated microphysical EOS. 
General relativistic hydrodynamic core-collapse si mulations 
using the ta bulated EOS were p erform ed by lOtt et all 
(2007a) and IDimmelmeier et al.l (I2007D. and magneto- 
hydrodynamic simulations by Cerda-Duran et all (|2007t ) 
using the passive magnetic-field approximation. These 
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three studies also included a simplified treatment of neu- 
trino transport. 

The rest of the paper is organized as follows. Section [2] 
presents a brief overview of the theoretical framework we 
use, namely the CFC equations and the GRMHD equa- 
tions in the 3-1-1 formalism. Our numerical approach is 
discussed in Sect.[3l Tests of the numerical code are pre- 
sented in Sect.U] including a magneto-rotational core col- 
lapse simulation, and the conclusions are given in Sect.[51 
Throughout the paper, we use a spacelike metric signature 
(— , -|-, and units where c = G = 1. We absorb the 
factor appearing in the MHD equations in the defi- 

nition of the magnetic field B^, i.e. the units of the magnetic 
field are ^/in Gauss. Greek indices run from to 3, Latin 
indices from 1 to 3, and we adopt the standard Einstein 
summation convention. 



2. Physical approach 

We adopt the 3 4- 1 formalism of general relativity 
(jLichnerowicd [194^ to foliate the spacetime into spacelike 
hypersurfaces. In this approach, the line element reads 



ds^ = -a^ de -f {dx' + dt){dx' + (3^ dt), 



(1) 



where a is the lapse function, /?* is the shift vector, and 7^- 
is the spatial three-metric induced in each hypersurface. 
Using the projection operator _L(^ and the unit four-vector 
normal to each hypersurface, it is possible to build the 
quantities 



E 



a 1 



1-t n''T^, = --{To.,-T,jl3^), 



(2) 
(3) 
(4) 



which represent the total energy, the momenta, and the 
spatial components of the energy-momentum tensor T^,y, 
respectively. 

To solve the gravitational field equations we choose the 
ADM gauge in which the three-metric can be decomposed 
as 7ij = 4''^^ij + hjj^ , where (j) is the conformal factor, 
7ij is the flat three-metric, and hj^ is the transverse and 
traceless part of the three-metric. We note that this gauge 
choice implies the maximal slicing condition where the trace 
K of the extrinsic curvature tensor Kij vanishes. 



2.1. The CFC approximation 

In our work, Einstein's field equations are formulated and 
solved using the conformally flat condition (CFC hereafter), 
introduced by Jsen bcrg (1978) and first used in a dynami- 
cal context bv I Wilson et al.l ([1996). In this approximation, 
the three-metric in the ADM gauge is assumed to be con- 
formally flat, 7ij = </'^7ij. We note that this approxima- 
tion can also be realized for other gauge choices such as 
the quasi-isotropic gauge or the Dirac gauge, both supple- 
mented by the maximal slicing condition. Under the CFC 
assumption, the gravitational field equations can be written 



as a system of five nonlinear elliptic equations. 



167r 



Idtt 



(5) 
(6) 



A/3' = 16™(/.45.^20iOif'^V, - iv'Vfc/3^ (7) 



where A and V are the Laplace and nabla operators asso- 
ciated with the flat three- metric, and S = ^/^^ Sij. 

2.2. General relativistic magnetohydrodynamics 

The energy-momentum tensor of a magnetized perfect fluid 
can be written as the sum of the fluid part and the electro- 
magnetic field part. In the so-called ideal MHD limit (where 
the fluid is a perfect conductor of infinite conductivity), the 
latter can be expressed solely in terms of the magnetic field 
h^^ measured by a comoving observer. In this case, the total 
energy-momentum tensor is given by 



(8) 



where p is the rest-mass density, h — 1 + e + P/ p the rela- 
tivistic enthalpy, e the specific internal energy, P the fluid 
pressure, the four- velocity of the fluid, and &^ = h^h^. 
We define the magnetic pressure Pmag = 6^/2 and the spe- 
cific magnetic energy e^ag = ^^/(2p), whose effect on the 
dynamics is similar to that of the fluid pressure and the 
specific internal energy of the fluid, respectively. 

For an Eulerian observer, = n^, and in the ideal 
MHD limit, the temporal component of the electric field 
vanishes, E^^ = where eijk is the permu- 

tation tensor and B'^ is the magnetic field. In this case. 
Maxwell's equations reduce to the divergence-free condi- 
tion and the induction equation for the magnetic field. 



0, 



dB* 



dt 



Vj{v*'B*^ -v*^B* 



(9) 



with B** = v^i?' and w* ' = av^ — /3% where is the fluid 
three-velocity as measured by the Eulerian observer. The 
ratio of the determinants of the three-metric and the flat 
three-metric is given by 7 = 7/7. 

The evolution of a magnetized fluid is determined by the 
conservation law of the energy-momentum, V^T'^'^ = 0, 
and by the continuity equation, V^J'^ = 0, for the rest- 
m ass current J'^ = pu '^- Following the procedure described 
bv I Anton et al.l (|2006f l. the conserved quantities are chosen 
in a way similar to the p urely hydrodynamic case presented 
bv lBanvuls et al.l (|l997t) : 



D = pW, 

5, = {ph + b^)W'^v^-ab,b°, 



{ph + b^)W' 



P 



a\b° 



(10) 
(11) 
(12) 



where W = au^ is the Lorentz factor. With this choice, 
the system of conservation equations for the fluid, and the 
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induction equation for the magnetic field can be cast as a 
first-order, flux-conservative, hyperbolic system. 



1 



(13) 



with the state vector, the flux vector, and the source vector 
given by 



U 

S = 



[D,Sj,T,B% 



0, iT^'^?^,a 
2 dx^ 



a- 



W 



W ' 



rptlO 9 hi a ^ rpfj.iy pO 



,0* 



(14) 

(15) 
(16) 



where i5* is the Kronecker delta, and r^_^ are the ChristofFel 
symbols associated with the four-metric. We note that the 
above definitions contain components of the magnetic field 
measured by both a comoving observer and an Eulerian 
observer. The two are related by 



WB^Vi 



V = 



B' + ab°i 
W 



(17) 



The hyperbolic structure of Eq. (I13|) and the associated 
spectral decomposition (into eigenvalues and eigenvectors) 
of the flux-vector Jacobians are given in lAnton et al.l 
(|2006( l. This information is required to numerically solve 
the system of equations using the class of high-resolution 
shock-capturing schemes that we have implemented in our 
code. 



2.3. Equation of state 

The new numerical code can handle a variety of equations 
of state including a polytropic EOS, an ideal gas EOS, a 
hybrid EOS, and a tabulated microphysical EOS. 



2.3.1. Hybrid EOS 

The hybrid EOS (jJanka et al.l I1993D is a simplified ana- 
lytical equation of state used in core collapse simulations. 
The pressure is given by a polytropic part, Pp = Kp'^, 
with K = 4.897 x 10^^ (in cgs units), plus a thermal part, 
^th = peth(7th — 1)) where the specific thermal energy, 
eth = e — Cp, and 7th = 1.5. The thermal contribution 
takes into account the increase of the thermal energy due 
to shock heating. When p exceeds the nuclear saturation 
density, Pnuc = 2.0 x lO^'* g cm~^, the value of 7 is raised 
to 72 = 2.5, and K is adjusted accordingly to guarantee 
the continuity of P and e. Due to this stiffening of the EOS 
the core undergoes a so-called pressure-supported bounce. 
More details about the hy brid EOS can be found, e.g. in 
iDimmelmeier et al.l (|2002a[ ). 

2.3.2. Microphysical EOS 

We further emplo y the tabu la ted n on-zero temperature 
nuclear EO S bv Shen et al.l (|l998l ) in the variant of 
iMarek et al.l (|2005t ). which includes baryonic, electronic, 
and photonic pressure components. It specifies the fluid 



pressure P (and additional thermodynamic quantities) as 
a function of p, the temperature T, and the electron frac- 
tion Yg. Whenever it is necessary in the code to compute 
the pressure as a function of the specific internal energy e 
instead of the temperature T, we iterate the corresponding 
value of T with a Newton-Raphson scheme. 



3. IMumerical methods 

Since our new numeri cal code is based on a previ ous purely 
hydrodynamic code (jPimmelmeier et al.ll2002"al lbl) and on 
its extension to the passive magnetic-field approxim ation 
(|Cerda-Duran fc Fond[2007t ICerda-Duran et al.ll2007f ). we 
describe here in detail only those numerical techniques that 
represent improvements over their predecessors, and pro- 
vide only concise information about the numerical schemes 
already described and tested elsewhere. 

The code solves the coupled time evolution of the equa- 
tions governing the dynamics of the spacetime, the fluid, 
and the magnetic field in general relativity. The equations 
are implemented in the code using spherical polar coordi- 
nates {t,r,6^ip}. We assume axisymmetry and equatorial 
plane symmetry. 



3.1. Metric solver 

The CFC metric equations, Eqs. ([5]|6|), are five nonlinear 
elliptic coupled Poisson-like equations, which can be writ- 
ten in compact form as Au{x) = f{x;u{x)), where u = 
u'' = ((/), a4>, (3^), and f — f'' is the source vector. These five 
scalar equations are coupled via the source vector, which de- 
pends on the components of u. We use a fix-point iteration 
scheme in combination with a linear Po isson solver to solve 
these equa tions (for further de tails see ICerda-Duran et al.l 
(|2005f ) and IDimmelmeier et all (2002a}). 

Since the CFC equations are written in terms of the 
energy-momentum tensor, the contribution of the energy- 
momentum of the magnetic field is automatically incorpo- 
rated in the computations. The main difference with re- 
spect to the non-magnetized case arises from the fact that 
the magnetic-field contribution does not necessarily have 
compact support. However, the magnetic field far from the 
fluid should decay at least as a dipole ^/r^), i.e. its con- 
tribution can be computed correctly by integrating the CFC 
equations in a sufficiently large volume. We checked that in 
all cases considered here the contribution of the outer mag- 
netic field to the energy-momentum tensor is too small to 
affect the CFC metric. Hence, we consider the contribution 
of the magnetic field only up to a radius ~ 20% larger than 
the radial extent of the fluid. The contribution of the inner 
magnetic field is, however, in some case sufficiently large 
to modify the metric (e.g. for the magnetized neutron-star 
equilibria), and therefore cannot be neglected in the CFC 
equations. 

3.2. Riemann solver 

For the evolution of the matter fields we utilize a HRSC 
scheme to integrate the subset of equations in the system 
of Eq. p3p that corresponds to the hydrodynamic variables 
(D, Si,T). HRSC schemes ensure the numerical conserva- 
tion of physically conserved quantities and a correct treat- 
ment of discontinuities such as shocks (see e.g. iFontI [20031 
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for a review and references therein). We implemented var- 
ious cell-reconstruction procedures that are accurate to ei- 
ther second-order or third-order in space, namely minmod, 
MC, and PHM (see iTorol [19991 for definitions). The time 
update of the state vector U relies on the method of lines 
in combination with a second-order accurate Runge-Kutta 
scheme. The numerical fluxes at cell interfaces are obtained 
using either the HLL single-state s olver of Harten et al 



1983 ) or the symmetric scheme of iKurganov fc Tadmoi 



20001 ) (KT hereafter). Both solvers yield results with an 



accuracy comparable to Riemann solvers exploiting the 
full characteristic information, as demonstrated for hy - 
drodynamic special relativistic (|Lucas-Serrano et al.ll2004 ) 
and general relativist ic flows in dynamical spacetimes 
(jShibata fc F ontl l2005D. Tests of both solvers in GRMHD 
were reported by I Anton "eTaLl (|2006h . 



3.3. Constrained transport scheme 

The evolution of the magnetic field needs to be performed 
differently from the rest of the conservation equations be- 
cause the physical meaning of the corresponding conserva- 
tion equation is different. Although the induction equation 
can be written in a ffux-conservative form, a supplementary 
condition for the magnetic field (the divergence constraint, 
or the conservation of the magnetic flux) has to be fulfilled 
during the whole evolution. A mong the numerical schemes 
that satisfy this condition (see lTothI i 



2000I. for a review), the 



constrained transport (CT) scheme (|Evans fc Hawlevlll988l ) 
was proven to be adequate for performing accurate simula- 
tions of magnetized flows. Our particular implementation 
of the CT scheme is adapted to the spherical polar coor- 
dinates used in the code, and uses cell interface-centered 
poloidal and (because of the assumption of axisymme- 
try) cell-cent ered toroidal magnetic-fiel d components (see 
Sect. 3.2.1 in ICerda-Duran fc Fontl[2007. for details). The 
induction equation is discretized in the same way as for the 
fluid equations. 

CT schemes preserve the magnetic flux during the evo- 
lution of a magnetized flow, but do not impose the di- 
vergence constraint on the initial magnetic field. Hence, 
one also has to provide initial data that fulfill this con- 
straint in order for the CT method to work properly. 
This can be ensured by computing the staggered magnetic 
field from the vector pote ntial (see Eqs. (28) and (29) in 
ICerda-Duran fc Fontl[2007l . for details). 

Finally, one has to consider the computation of cell- 
centered values of the (poloidal) magnetic field, which are 
required in the source terms and for the reconstruction of 
the magnetic field tangential to the cell interfaces, that 
enter the evaluation of the num erical flux. Here, w e de- 
part from the scheme described bv lAnton et al.l (|2006t) . who 
computed cell-centered magnetic fleld components assum- 
ing that the corresponding magnetic flux at the cell center 
is given by the average of the magnetic flux at the cell in- 
terfaces. Using this prescription, the cell-centered magnetic 
pressure differs from that at the interface, even in the case 
of a homogeneous magnetic field. Instead we use 



B* 



B* 



cos 6 



cos 9 



3+i 



COS 6, 



* J -2 

sm I 



B* 



sin 6. 



sin 1 1 



« 3- 



sm( 



(18) 



for the cell-centered magnetic-field components. This pre- 
scription guarantees that for a homogeneous field parallel 
to the rotation axis, the magnetic pressure is equal at the 
cell center and the cell interface. It also increases the stabil- 
ity of the code for highly magnetized fiows, especially near 
MHD equilibria, and is critical for the success of some of 
the tests presented here. 

3.4. Recovery of primitive variables 

In relativistic hydrodynamics, in contrast to the Newtonian 
case, there exists no explicit expression for the primitive 
variables (p, , e) in terms of the conserved ones {D,Si,T). 
Hence, a recovery procedure is required whereby the prim- 
itive variables are obtained from the conserved ones by in- 
verting the nonlinear system given by Eqs. (|10l- ll2p with an 
efficient numerical algorit hm. In most of the recovery algo- 
rithms (jNoble et al.ll2006t ). one first introduces some scalar 
quantities, and then solves the resulting simplified system 
of equations before recover i ng th e primitives. 

Following lAnton et all (|2006t) . our recovery procedure 
is based on the two scalar quantities (note that the first of 
these is the conserved energy) 



(19) 



T = phW"^ -P + b^{W^ - 1/2) - a^{lPf 

+a\b°)^ [-2phW^ + b\l - 2W^) + a2(60)2] . (20) 

If one defines z = phW'^ and makes use of the expression 
B ■ S = phWalP to eliminate &° from these equations, the 
resulting expressions 



2z + B2 



{B ■ S)' 



-{z + B^f = 0, (21) 



D 



B' 



{B-S) 

2z2 



p 



7 



0, 



(22) 



depend only on conserved quantities, on the metric, and 
on the set of unknowns {P, z, W}, respectively. The system 
formed by Eqs. (PT] - [^ and the EOS can then be solved to 
obtain {P, z, W}. From these three quantities, the primitive 
variables can be easily computed as 



P 



D 

-/'^Sj + (B ■ S)B'/z 

z - DW - PW^ 
DW ■ 



(23) 
(24) 
(25) 



The numerical procedure to solve the system of Eqs. ((2T] - 
[22]) therefore depends on the EOS (see next two subsec- 
tions). 

3.4.1. Barotropic fluid 

In a barotropic ffuid, the pressure depends only on the den- 
sity, i.e.P{p). In many astrophysical situations (e.g. cold 
neutron stars) as well as in many standard tests of numer- 
ical codes, the fluid is assumed to be barotropic. The most 
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commonly used barotropic EOS is the polytropic EOS, 
P = Kp^, where K is the polytropic constant and F is 
the adiabatic index. 

For a barotropic EOS, the enthalpy is a function of the 
density only, i.e.h{p), and thus z = Dh{p)W. Using this 
fact and Eq. ((23)) . it is possible to eliminate the unknowns 
P and z from Egs. (PT] - [^ . The Lorentz factor W then 
remains to be computed numerically by solving one of these 
equations. 

Following lAntonI ()2007[ ). we solve Eq. ^ for W by 
means of the bisection method, and then recover P and 
z using the EOS. This method is extremely robust and al- 
ways leads to a solution for W, provided that it lies between 
the initial lower and upper guess value. 

3.4.2. Baroclinic fluid 

This is the most common form of the EOS in hydrodynamic 
simulations because it takes into account temperature ef- 
fects. We implemented several baroclinic EOS in our nu- 
merical code, namely t he ideal gas EOS, P = pe(F — 1), the 
analytic hybrid EOS (jjanka et al.lll993h . and a tabulated 
microphysical EOS. 

Irrespective of the baroclinic EOS used, it can always be 
expressed in the form P{p, e,Yi). Since the composition Yi 
(the index i runs over all relevant species) is known directly 
from the hydrodynamics, the dependence of the EOS on the 
composition does not affect the recovery procedure. The 
following discussion therefore can be restricted to an EOS 
of the form P{p, e). 

For a baroclinic fluid, the system formed by Eqs. pTl - 
and the EOS expressed as 

P-P(p,e)-0 (26) 

must be solved numerically. In general, it is not possible 
to use the EOS to analytically remove the dependence of 
Eqs. (|2T] -[22 | on P, since e depends on the pressure itself due 
to Eq. (P5)l . However, for some analytic EOS, e.g. an ideal 
gas, Eqs. (^5]) and (^5]) can be used to express the pressure 
as a function of z and W only. This allows one to eliminate 
P from Eqs. (|21l - [^^ . reducing the system to be solved to 
two equations with the unknowns z and W. Since the nu- 
merical method should not rely on any assumption about 
the EOS, we do not exploit this simplification. Instead we 
solve the system of Eqs. (dU [221 IHl) by means of a Newton- 
Raphson scheme, which converges rapidly provided the ini- 
tial guess is sufficiently good (see below). We consider the 
Newton-Raphson iteration to be converged when the rela- 
tive error of the variables is less than a certain tolerance 
(typically IQ-^^). 

Microphysical tabulated EOS. Some of the equations of state 
available from nuclear physics that are used in astrophysics 
are not provided in terms of the specific internal energy. In 
general, the EOS depends on the composition of the fluid 
(usually the electron fraction Yg) and on the temperature 
T instead of e. Hence, one has to deal with tabulated EOS 
of the form 



p = P(p,T,ye), 

e = e{p,T,Y,). 



(27) 
(28) 



A first approach to handle such an EOS is to obtain ef- 
fectively P = P(p, e. Ye) by computing the value of T that 



satisfies Eq. ([28)1 for a given value of e. The procedure de- 
scribed above for an EOS of the form P{p, e) can then be 
applied, since Y^, is known directly from the evolution. This 
approach was succes sfully used in the hydrodynamic simu- 
lations presented bv lOtt et al.l ([2007bl ). In the magnetized 
case, however, we find this approach to be problematic for 
strong magnetic fields {Pmag/P > !)■ The solver is able to 
recover the exact value in that regime only if the initial 
guess is very close to the solution, which renders the code 
unstable. To avoid this problem we add the equation 



e-e(p,r,ye) -0, 



(29) 



to the Newton-Raphson system, and solve the extended sys- 
tem of Eqs. ([m [11 [Ml) for the unknowns z, VF, and T. This 
allows one to use directly the EOS as a function of T in- 
stead of e. We find that this method is very stable and has 
a much larger radius of convergence than the first approach 
(see Sect. liTT]) . 

"Safe" guess values. When we use the values of the previous 
time step as an initial guess for the Newton-Raphson itera- 
tion at a given time step, the solver usually converges within 
a few iterations. However, sometimes the guess values are 
too far away from the solution, and the Newton-Raphson 
iteration fails. In such a case, we restart the iteration pro- 
cess using a "safe" set of guess values, which we choose to 
be upper limits to the unknowns {P, z, W} (or {T, z, W}). 
This choice leads to a rather robust recovery scheme as 
demonstrated by our test calculations (see Sect.[l]). 

To derive an upper limit for z, we define i5 as the angle 
between v and B, i.e. v B = \J v'^B'^ cos S. Using this angle, 
we have 



b' = 



El 



(l + [W^ - l)cos2 5) 



From the definition of r given in Eq. (|12p . one obtains 
~2 



P + D 



-(2- 



(30) 
(31) 

z as 

(32) 



2 

Since the last term in this equation is always negative or 
zero, an upper limit for z is given by 



z <T + P 



2 



(33) 



However, this upper limit cannot be computed directly from 
the conserved quantities, since the pressure is unknown. 
Hence, we first need to determine an upper limit for the 
pressure. If we assume that the pressure grows monotoni- 
cally with p and e, which is a reasonable assumption for the 
types of EOS that we use in the code, then we only need 
to derive upper limits for p and e, and hence 

P — Pmax — P(Pmax; ^max)- (^^) 

It is easy to find an upper limit for p, since W > I, 

P < Pmax = D. (35) 

In the case of the specific internal energy, we substitute 
Eq. into Eq. ([25]) and obtain 



1 

DW 



B 

T - — + P(l - iy2) _^ JJ{1 _ p^) 



-(2 



cos' 



'5) 



(36) 
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Using again the fact that > 1 we derive the upper hmit 



1 
D 



An upper Hmit for z is given by 



(37) 



(38) 



which coincides with z in the hmit of smaU velocities, W 
1. 

We were unable to compute an analytic upper limit for 
W; however it is easy to set an upper limit from physical 
considerations. In the core-collapse simulations in which we 
are interested, the Lorentz factor is not expected to exceed 
a value of 10. Nevertheless, we chose a much larger guess 
value for W, since the number of iterations until conver- 
gence is very insensitive to the precise value of the upper 
limit. Accordingly, the "safe" guess values that we use in 
the Newton-Raphson solver are 



P — P 

^ guess — max 5 



•^gucss — ^max i 
M^^uess = 10000, 



^gucss — ^i,P 



max; ^inax 



(39) 
(40) 
(41) 
(42) 



We note that Tguess is not an upper limit for the temper- 
ature in general, but the pressure value computed with 
and Pmax provides an upper limit for the pressure. 



^ guess 



i.e.P < P{p 



^(Tunss ) . 



max; J- guess 



3.5. Vacuum treatment 

The presence of vacuum regions is common in numeri- 
cal simulations dealing with astrophysical scenarios. These 
regions are usually avoided by imposing a numerical at- 
mosphere surrounding the object under study, i.e. a small 
floor value for the rest mass density, which allows one 
to use the same recovery procedure in regions filled with 
the numerical atmosphere and the fluid. A vacuum re- 
gion would cause the recovery procedure to fail both in 
the hydrodynamic and in the magneto-hydrodynamic case, 
as can be inferred from Eqs. (|24l[^5)) . The numerical at- 
mosphere approach is commonly used in hydrodynamic 
simulations (see e.g. iFont et aLll2002t iDimmelmeier et al.l 
2002al) as well as in GRMHD simulations (|Du ez et al."2005|; 
Shibata k Sekiguchi|[200l iGiacomazzo fc R czzoUa 2007). 

In the unmagnetized case, the floor value of the numeri- 
cal atmosphere is chosen such that it does not affect signif- 
icantly the dynamics of the system. This can be achieved 
by choosing the threshold value for the rest mass density 
to be a small fraction of the maximum density in the ini- 
tial model, typically pthr ~ 10~®Pmax- Every grid point 
with p < pthr is reset to the numerical atmosphere value, 
i.e. p = Patm and w* = 0, where the floor value for the rest 
mass density is Patm ~ 10~'^pthi - 

In the magnetized case additional problems arise. Since 
the transition to the numerical atmosphere usually results 
in a steep profile in p (which drops to the floor value) but 
not necessary in B (magnetic field lines can extend into the 
vacuum) , atmosphere regions can easily have large ratios of 
Pma.gl even if the fluid is weakly magnetized. This prob- 
lem increases as the floor value patm is reduced, and can lead 
to problems with the recovery of the primitive variables 



in the atmosphere. To avoid this problem, some authors 
(|Duez et al.ll2005t IShibata fc Sekiguchil [20051) do not allow 
magnetic fields in the numerical atmosphere by choosing 
magnetic fields confined to the fluid regions. Other authors 
(jGiacomazzo fc RezzoIIall2007tlShibata et al.ll2007D apply a 
floor to the hydrodynamic variables and allow the magnetic 
field to evolve freely. This approach works fine, if the ratio 
of Pmag/P does not exceed the critical value above which 
the recovery procedure fails. We estimate this critical value 
for our code in Sect. 14. II Consequently, a sufficiently low 
density atmosphere will show the correct dynamic behavior 
when the magnetic field strength in the atmosphere is lim- 
ited. On the other hand, if one wishes to simulate stronger 
magnetic fields, one must use a denser at niosphere that 
can e ven affect the dynamics of the system ([Shibata et al.l 
120011) . 

We also allow for a freely evolving magnetic field in the 
atmosphere, since we are then not restricted to any par- 
ticular magnetic field structure. To overcome the problem 
of the high magnetization Pmag/f in the atmosphere, we 
choose compromise values for patm depending on the prob- 
lem to be solved. Since there are cases (e.g. for the evolution 
of neutron stars in Sect. 14. 3i Pmag/P becomes as large as 
10"'^^) where the ratio Pmag/P exceeds the critical value for 
the recovery procedure (10® — 10®; see Sect. 14. ip even for 
reasonable values of patm , we use a fast atmosphere check- 
ing routine which avoids the recovery of the primitives in 
the respective zones. 

If we are able to mark a zone as being part of the atmo- 
sphere before the recovery of the primitives is performed, 
we can avoid the recovery because the values for the primi- 
tives in these zones are set to the floor value. Since D > p, 
if D < pthr holds, then p < pthr, and the zone is part of 
the atmosphere. Hence, we can use this condition to check 
whether a zone belongs to the atmosphere before perform- 
ing the recovery. We note that for atmosphere zones, whose 
velocities are set to zero at every time step, it is very un- 
likely that the value of the (unknown) Lorentz factor W at 
the next time step differs significantly from 1, i.e. p ^ D in 
these zones. Using this procedure, we can handle arbitrarily 
large magnetic fields in the atmosphere without imposing 
any limitation on the value of 

Patm ■ 

In the equilibrium models of magnetized neutron stars 
f Sect. 14. 3) ). we keep the magnetic field fixed in the atmo- 
sphere to its initial value. This is a reasonable choice since 
for an equilibrium configuration the outside magnetic field 
should not change during the evolution. The advantage of 
this approach is that the time step is not dominated by the 
atmosphere where the eigenvalues are close to the speed 
of light, but by the neutron star interior with eigenval- 
ues of the order of « 0.1. This results in a speed-up of 
about a factor of 10 in these computations. We note that 
this speed-up is only possible because in the CFC approxi- 
mation the metric evolution does not constrain the size of 
the time step. In codes based on hyperbolic formulations 
of the Einstein equations, the time step is always limited 
by the light-crossing time of the zones, i.e. this speed-up 
is impossible. All other existin g GRMHD codes with dy- 
namic spacetimes (iDuez et al.i 20 05: Shiba ta fc Sekiguchil 
120051 : [G lacomazzo fc Rezzollall20011) suffer from this limita- 
tion. 
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Table 1. Test of the recovery of the primitive variables. 
Varying the flow velocity v, the magnetic field B, and the 
angle between B and v in a, wide range, the recovery pro- 
cedure is tested for the equations of state, the densities p, 
the specific internal energies e, and the electron fractions 
Yg given in columns 2 to 5, respectively. 
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4. Code tests 

4.1. Recovery of the primitive variables 

The numerical method used for the recovery of the 
primitive variables {p,v^ , e, B^} from the conserved ones 
{D, Si, T, B^}, is tested by varying the flow velocity v, the 
magnetic field B, and the angle S between B and v over a 
wide range. Instead of varying v and B we vary the Lorentz 
factor — 1 in the interval [10~^, iC], and the magneti- 
zation Pmag/-P in the interval [10~^, 10^°], respectively. We 
choose values of p, e, and (for the tabulated EOS only) 
Ye that are typical for core collapse (Table [1]). Test cases 
PN and PE correspond to a polytropic EOS with F = 2, 
K = 1.455x10^ (cgs units), and P = 4/3, K = 4.897x10", 
respectively. For the hybrid EOS, the test cases HN and 
HE probe densities above and below nuclear matter den- 
sity, while the SHEN EOS cases test the typical conditions 
inside a proto-neutron star (SI) and the progenitor core (S2 
to S4). 

For a baroclinic EOS, we use the "safe" values given 
in Sect. 13. 4.^ as guess values for the Newton-Raphson it- 
eration, and choose a tolerance of 10~^^. Figure [T] shows 
the region in which the recovery scheme converges well, 
i.e. where the relative difference between the values of the 
recovered primitive variables and their exact values is less 
than 10"^". For all considered EOS, the parameter space of 
astrophysical interest is well covered. In test case 17, the re- 
gion of convergence is reduced substantially when cos S = 0. 
This test, that corresponds to an extreme ideal gas with 
e = 10^, was chosen to determine the maximum value of 
e which can be recovered and lies outside the parameter 
range of interest in core collapse. We also determine up- 
per limits to the magnetization Pmng/P in the low velocity 
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Fig. 1. Upper limits for the Lorentz factor {W — 1) and 
the magnetization Pniag/-P for which the relative difference 
between the values of the recovered primitive variables and 
their exact values is less than 10^^". The upper panel shows 
these limits for the polytropic and the hybrid EOS, the mid- 
dle panel for the ideal gas EOS, and the lower panel for the 
tabulated SHEN EOS, respectively. Thin lines correspond 
to the case cos^ 6 = 0, and thick lines to cos^ 5 = 1.0. The 
shaded region in the bottom panel shows the typical pa- 
rameter space encountered in core-collapse simulations. 

limit ranging from values of 10® to 10*. In the low magnetic 
field limit the maximum Lorentz factor that the recovery 
procedure can handle is 10^ to 10"^. If any of these lim- 
its is exceeded, the numerical scheme is unable to recover 
the primitive quantities within the required accuracy. The 
reason for the recovery failure in the three limiting cases 
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for e, W, and Pmag/P is that the contribution of internal 
energy, kinetic energy or magnetic energy, respectively, is 
dominant in the system and any other kind of energy has 
a very small contribution. In these cases large changes in 
the subdominant terms will produce small changes in the 
recovery equations and the system may therefore converge 
to a wrong solution within a given accuracy. If the toler- 
ance value is reduced, these limits can be extended. We note 
that in astrophysical situations involving baryonic matter, 
it is unlikely to encounter values oi e > 1, W > 10, or 
Pm&g/P > 100. Therefore, we consider that our recovery 
procedure is sufficiently robust for simulations of core col- 
lapse and involving compact objects. 

Using the "safe" guess values, the number of iterations 
needed for the Newton-Raphson solver to converge is rela- 
tively large: 50 — 70 for both the hybrid and the ideal gas 
EOS, and 50-200 for the tabulated EOS. However, during 
a numerical simulation, the "safe" guess is only used if the 
regular guess (the value from the previous time step) fails. 
If we use guess values that differ by only 10% from the exact 
ones, the Newton-Raphson converges more rapidly, within 
10 — 20 iterations for the hybrid and ideal gas EOS, and 
20 - 30 for the tabulated EOS. For the polytropic EOS it 
takes about 40 — 60 bisection steps to achieve the required 
tolerance. 

4.2. Spherical explosion 

Since the majority of existing (2D) MHD codes are written 
in cylindrical coordinates, a commonly performed test is the 
simulation of a cylindrical explosion. For relativistic MHD 
codes, such a setup was proposed by iKomissarovl (| 19991) . 
which was also used by o ther authors ( Del Zanna et al.l 
120031 : iLeismann et al.l [20051 ). However, when using a code 
based on spherical c oordinates the mo st natural choice is a 
spherical explosion. iKossl et al.l ()1990D performed this test 
with a Newtonian MHD code, but to the best of our knowl- 
edge no spherical explosions test has been performed in 
relativistic MHD. Therefore, we consider here a spherical 
explosion test for which the initial jum p conditions are iden - 
tical to those of the cyhndrical test of IKomissarovl ()1999D . 

Our test setup consists of an initial radially symmetric 
explosion zone (r < 1) surrounded by a highly magnetized 
ambient gas for r > 1. In the outer part of the explosion 
region (0.8 < r < 1.0), we set the state variables decline 
exponentially to the values of the ambient medium (Table 
[2]). The velocity is initially zero everywhere, and the mag- 
netic field is homogeneous and parallel to the symmetry 
axis. The background spacetime is assumed to be flat. The 
initial data are evolved using an ideal gas EOS with an 
adiabatic index T — 4/3. The computational grid is evenly 
spaced in radius and angle, and extends in the radial direc- 
tion up to a maximum radius of r = 6.0. We perform the 
test with two (r, 6) resolutions (80 x 20, and 160 x 40) for 
all reconstruction schemes and flux formulae. 

Table [5] shows the eigenvalue structure of the initial 
setup. Since the explosion region is weakly magnetized 
{Pmag/P = 5 X 10"'^), the dominant wave in this region 
is the fast magnetosonic wave, which propagates at a speed 
I A/ 1 close to that of the corresponding hydrodynamic wave. 
All other wave speeds are close to zero, i.e. the inner re- 
gion will expand with the velocity |A/|. The ambient gas 
is highly magnetized, and the full wave structure is signif- 
icant there. The only waves fast enough to travel ahead of 



the explosion shock are the fast magnetosonic wave, which 
propagates with an almost angular-independent radial ve- 
locity close to the speed of light, and the Alfven wave, whose 
radial velocity is also close to c along the symmetry axis, 
but varies as cos 61. 

Figure [2] shows a snapshot at the end of the simulation 
{t — 4). Both the Lorentz factor and the pressure distri- 
bution clearly show the wave structure mentioned above. 
The spherical fast magnetosonic wave is located at r « 5, 
and the trailing strong shock, which is deformed due to the 
magnetic field, consists of a mixture of the bulk expansion 
of the inner region and an Alfven wave propagating faster 
along the axis. Even further inwards, a rarefaction wave 
is visible, which is almost spherically symmetric since the 
magnetization in this region is rather low. 

The corresponding radial profiles of P and W both 
along the equator (upper panels) and the axis (lower pan- 
els) are displayed in Fig.|3| for various numerical methods. 
These plots are qualitatively similar to those of the cylindr i- 
cal explosion test (see e.g. Fig. B.4 in ILeismann et al.ll2005D . 
All numerical schemes exhibit first order convergence as ex- 
pected for flows with shocks. The MC and PHM schemes 
yield very similar results, while the minmod scheme gives 
slightly smaller values. No significant differences are found 
between the results obtained using the HLLE and KT fiux 
formulae. 

4.3. Magnetized neutron stars 

The previous two tests demonstrate the ability of the 
code to handle extreme situations such as high magneti- 
zation, large Lorentz factors, and strong shocks. In this 
section we show its correct behavior in curved space- 
times, particularly in dynamic ones. An astrophysical sce- 
nario that can be used for this assessment is the evo- 
lution of equilibrium neutron stars, a test which is fre- 
quently used fOTgene^^ 

(Shibata'1999; 'Font et al."2002'; 'DimmclmeiCTeLalJl2002i 
Ducz ct al. 2003- Cerda-Duran_ct_al, 20Q^ as well as 
for GRMHD codes (jGiacomazzo fc Rezzollal 1200 Tt) . Of all 
presently existing codes capable of solving th e GRMHD 
equat i ons coupled to a dynami c spacetime (iDuez et al 



2005'; 'Shibata & SekiguchJ 120051: [Giacomazzo & Rezzolla 



2007: .Anderson et al. 2008D . this demanding test, involving 
all aspects of the code and in particular the correct cou- 
pling between metric and MHD equations, has only be en 
performed by the code of I Giacomazzo fc Rezzollal (I2007D . 

As initial models for the magnetized neutron star test, 
w e use the relat i vistic self-consistent equilibrium models 
of iBocquet et al.l (|l995l ). where all effects of the magnetic 
field (Lorentz force, spacetime curvature generated by the 
magnetic contribution to the energy-momentum tensor) 
are taken into account. The equilibrium models are com- 
puted using the LORENE library Q. We construct non- 
rotating polytropic equilibrium models with F = 2 and 
K = 1.455 X 10^ (cgs units). The central enthalpy is chosen 
to be In/ic = 0.228, and the magnetic field i s that of a per- 
fect c onductor with the current density of iBocauet et al.l 
(I1995D and vacuum outside. By increasing the value of the 
central current density jo from to 5 x lO^^Am"^, we 
compute a sequence of equilibrium models with a magnetic 
field ranging from zero to 1.8 x 10^^ \/47rGauss (Table [3]). 



http:/ /www. lorene.obspm.fr/ 
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Table 2. Spherical explosion test. Initial values for the pressure P, the density p, the magnetic field and the 
magnetization Pmng/P are given in columns 2 to 5 for the explosion region (r < 1) and in the ambient region (r > 1), 
respectively. The eigenvalues, namely the speeds of the fast magnetosonic wave A/±, the Alfven wave Xa±, the slow 
magnetosonic wave Xs± , and of the entropy wave Ae in the radial direction at the initial time are given in columns 6 to 
9, both for the equator (left value) and the pole (right value). In addition, we provide in the last column the value of the 
eigenvalue (sound speed) A± of the corresponding non-magnetized case. 
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Fig. 2. Snapshot of the spherical explosion test at t = 4. The panels show the logarithm of the pressure (left) and 
Lorentz factor, and the magnetic field lines (right). The simulation was performed with 160 x 40 zones using PHM cell 
reconstruction and the KT flux formula. 



Table 3. Initial models of magnetized neutron stars. From left to right, the columns give the central current density jo, 
the equatorial radius r^, the ratio of polar to equatorial radius rp/r^, the ratio of magnetic to thermal pressure Pmag/P 
at the center of the star, the central magnetic field \B\c, and the ADM mass Madm of each model, respectively. 
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The magnetic field topology is shown in Fig. 2] for a rep- 
resentative model (MNS3). It is purely poloidal with field 
lines crossing the surface of the neutron star (thick dash- 
dotted line). At sufficiently large distances from the star, 
the magnetic field has a dipole topology. 

First, we perform simulations in the Cowling approx- 
imation, where the spacetime is kept fixed. We stop the 
evolution after 5 ms which corresponds to 52tdyn7 where 
^dyn = ^/rf/M is the characteristic dynamic time-scale of 
the system. Using the Cowling approximation, allows us 
to test the behavior of our MHD scheme without including 
yet the coupled evolution of the spacetime itself. The space- 
time fields are computed using the CFC equations in the 
first time step, and their values are kept fixed afterwards. 
To carry out convergence tests, we performed computations 
with models MNSO and MNS3 on equidistant grids (n^ x ng) 
with 80 X 10, 160 x 20, and 320 x 40 zones, respectively The 
other two models, MNSl and MNS2, were simulated only 



with 160 X 20 zones. We use the PHM reconstruction scheme 
and the KT flux formula in all computations reported in 
this section. The neutron star is surrounded by an atmo- 
sphere as described in Sect. 13. 51 with a threshold value of 
Pthr = lO-'^/Omax, and a floor value patm = lO'^Pmax- In 
the highly magnetized models MNS2 and MNS3 the value 
of Pmag/P is close to the critical value for the recovery pro- 
cedure in the outermost zone of the neutron star. In these 
models, we raise the threshold value to pthr = lO-^Pmax 
keeping the same floor value. 

Three of the panels of Fig. [5] show the evolution of the 
central density with time. Due to numerical truncation er- 
rors in the remapping of the equilibrium model from the 
spectral grid used by LORENE to our finite-difference grid, 
some small amplitude perturbations are triggered, which 
excite the normal modes of pulsation of the star. This 
causes the periodic oscillations of the central density. The 
neutron star remains in equilibrium throughout its evolu- 
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Fig. 3. Results for the spherical test explosion at < = 4. The panels show radial profiles of the fluid pressure P (left panels) 
and the Lorentz factor W (right panels) along the equator (upper panels) and the polar axis (lower panels), respectively. 
The lines styles differentiate between the reconstruction schemes: minmod (dotted), MC (dashed), and PHM (solid). 
Results for two different flux formulae are shown: HLLE (thin lines) and KT (thick lines). Note that the results obtained 
with these two flux formulae are often so similar that they cannot be distinguished. The grid resolution used is 160 x 40 
zones. 



tion, only a small drift with time is visible in the central 
density evolution. As we increase the grid resolution, this 
drift tends to zero. Comparing the non-magnetized model 
(MNSO) with the magnetized models (MNSl, MNS2, and 
MNS3), the drift, although small (< 0.2% in the 160 x 20 
models) , is larger in the magnetized models (lower left panel 
of Fig. [5]). We find that the drift is very sensitive to the value 
of pthr in model MNSO, if pthr > 10^^ Pmax (for smaller 
values, there is no influence). We suppose that, for denser 
atmospheres there, is a coupling between the star and the 
atmosphere that allows a transfer of mass and momentum 
from the interior to the atmosph ere ( see the related dis- 
cussio n in lStergioulas et all (|2004l ). and lDimmelmeier et al.l 
(|2006D ). This has two consequences for the evolution: first, 
the oscillations are damped more quickly, and second, the 
slope of the drift changes, even becoming negative. In the 
magnetized case, even if the magnetic field is weak, we have 
an extra coupling of the interior with the atmosphere due 
to the magnetic field lines leaving the star's surface. This 
causes an additional very small transfer of mass and mo- 
mentum from the atmosphere to the neutron star, which 
increases the drift in the evolution (see Fig. [5]). 

The convergence tests show that the order of conver- 
gence is 2.13 and 1.56 for model MNSO and MNS3, respec- 
tively. This global order of convergence is consistent with 
the second-order accuracy of our numerical TVD scheme. 



which reduces to first order at local extrema such as the 
center of the star and its surface. 

We also compute the Fourier transform of the central 
density evolution to obtain the mode frequencies of the neu- 
tron star pulsations (lower right panel of Fig. O. We find 
the fundamental mode frequency at about / = 2.7kHz, 
and subsequent harmonics at 4.6, 6.4, 8.2, 10.0, 11.8, and 
13.8 kHz, respectively. Since the energy of the magnetic field 
is small compared with the potential energy of the star, the 
influence of the magnetic field on the mode frequency is 
small. We find no frequency difference between the neutron 
star models within the frequency resolution (~ 0.5 kHz). 
We further observe that the quality of the spectrum de- 
teriorates at higher frequencies for models with stronger 
magnetic fields . We suspect that this degradation is an ar- 
tifact due to the stronger coupling of the interior with the 
atmosphere in the magnetized case. 

The second part of the test consists of the evolution 
of the same neutron star equilibrium models in a dy- 
namic spacetime. For reasons of computational efficiency, 
the CFC equations are computed only every 100th time 
ste p, the metric being interpola ted in-between as described 
by iDimmelmeier et al.l (jlOOli). The resuhs (Fig.© are 
qualitatively the same as those of the Cowling case dis- 
cussed before. The dynamic spacetime causes larger per- 
turbations in the central density evolution, which now also 
exhibits a larger drift with time (< 10% for the 160 x 20 
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Fig. 5. Evolution of equilibrium configurations of neutron stars in the Cowling approximation (fixed spacetime). The 
upper panels show the central density normalized to its initial value for the non-magnetized neutron star MNSO (left) 
and the magnetized model MNS3 (right). The results are displayed for three different grid resolutions (n^ x ng): 80 x 10 
(dotted), 160 x 20 (dashed), and 320 x 40 (solid). The lower left panel shows the evolution of p/ Pc o for all computed 
models for a grid resolution of 160 x 20: MNSO (sohd), MNSl (dashed), MNS2 (dotted) and MNS3 (dash-dotted). The 
lower right panel gives the corresponding Fourier transforms. 
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Fig. 4. Magnetic field lines structure of model MNS3. The 
thin solid lines represent the magnetic field lines while the 
thick dashed lines are rest mass density isocontours for 1, 
3, 5 and 7 x 10^"* g cm~^. Moreover, the thick solid line 
represents the surface of the star, and the thin dash-dotted 
line marks the boundary of the numerical grid. 



models). Similar drifts were already observed in fully cou- 



pled simulations of non-m agnetized (Font et al.l l2002f ) and 
magnetized models (jGiaco niazzo k Rczzolla 200^ In both 
models, MNSO and MNS3, the drift reduces with increasing 
resolution, and the order of convergence is 3.1 and 2.5 re- 
spectively. The convergence order is higher than expected 
(second order). We suspect that this is because the 80 x 10 
zone model is poorly resolved, i.e. the accuracy tends to 
grow faster than the order of convergence when doubling 
the resolution. Regarding the comparison between magne- 
tized and non-magnetized models (lower left panel of Fig^ , 
we observe larger drifts in the magnetized case due to the 
stronger coupling with the atmosphere. 

The Fourier transform of the central density for the 
160 X 20 models with dynamic spacetime evolution (lower 
right panel of Fig. [6]) gives a fundamental frequency of 
/ = 1.4 kHz, and higher harmonics at 4.0, 6.0, 7.8, 9.8, 
11.6 and 13.7kHz, respectively. We find no dependence on 
the amount of magnetization within the frequency resolu- 
tion. A similar resul t was obtained in the simulations of 
iMontero et al.l (|2007D regarding pulsating and magnetized 
thick accretion tori around Schwarzschild and Kerr black 
holes. This is unsurprising since the normal modes of a 
star are basically sound waves propagating in the radial 
direction, and the speed of sound is hardly altered by the 
magnetization of the investigated models. However, in the 
magnetized case, new modes can appear due to the richer 
eigenvalue structure of the GRMHD equations. In particu- 
lar, it is important to note that Alfvcn modes can be excited 
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in the star. For the magnetic field strengths present in our 
models, these mode frequencies lie below 100 Hz, i.e. much 
longer simulations are required to be able to see them in the 
spectrum. A deeper study of the Alfve n modes performed 
with o ur numerical code can be found in lCerda-Duran et al.l 
((2008h . 

If we compare the frequencies with those in the Cowling 
approximation, we observe that the Cowling approximation 
tends to overestimate the frequency of the modes (by al- 
most a factor 2 for the fundamental mode). The higher the 
order of the harmonics, the smaller is the overestimation, a 
trend that was observed before in numerical sim ulations of 
purely hydrodynamic models (|Font et al.ll2OO20 . The rea- 
son for this behavior is that perturbations on time scales 
smaller than the typical time scale of variations in the grav- 
itational field (which is roughly idyn) behave similarly as in 
a fixed spacetime. The frequency corresponding to the dy- 
namic time scale is /dyn — 10.4 kHz. Therefore, modes of 
frequency higher than /dyn will be unaffected if the com- 
putation is carried out in the Cowling approximation. This 
agrees with our mode computations. 



4.4. Core collapse 

The final test of our numerical code concerns simulations of 
magneto-rotational core-collapse. We note that these sim- 
ulations are not intended to be of astrophysical relevance, 
since the treatment of neutrinos in the code is still too 
poor for a study of the supernova explosion mechanism. 
Nevertheless, the tests allow us to validate the code in a 
fully dynamic context including strong magnetic fields, re- 
alistic stellar progenitors, and a microphysical EOS. To the 
best of our knowledge such demanding simulations have not 



yet been performed, which highlights the unique potential 
of our new numerical code for the study of relativistic stel- 
lar core collapse. 

As an initial model, we employ the inner part of the 
ir on core of the s olar-m etallicity 20 Mq progenitor model 
of' Wooslev et al.l (2002) . To this spherically symmetric and 
non-magnetized model, we add a rotation profile and a 
poloidal magnetic field. The rotation law for the specific 
angular momentum is given by j = A^{flc — ft), where 
A = 5 X km and fl is the angular velocity, which has 
a value ilc = 4.035 s~^ at center. The magnetic field is 
generated by a circular current loop of radius 400 km. This 
corres ponds to model s20AlB5-D3 in ICerda-Duran et al] 
(|2007) where a more detailed description can be found. 
We perform simulations for two different initial magnetic 
field strengths, namely for the weakly magnetized model 
S20A1B5-D3M10 with a central magnetic field of \B\c = 
10^° -s/Itt Gauss, and for the strongly magnetized model 
S20A1B5-D3M12 with |S|c = lO^^ Gauss. The models 
are evolved with the tabulated EOS of IShen et al] (Il998l) 
and an approximate dele ptonization scheme ( Liebendorfeil 
2005*) a s described by Dim melmeier et al] ( 20071 ). and 



Cerda-O Tuin et all ([2P07). We compare the evolutions of 
these two models with that of the co r respo nding model 
S20A1B5-D3M0 of ICerda-Duran et all (|2007[ ). which was 
evolved with the passive field approximation. Since the ef- 
fect of the magnetic field on the collapse dynamics is ne- 
glected for model s20AlB5-D3M0, its evolution should be 
similar to that of our weakly magnetized model s20AlB5- 
D3M10. The comparison with the passive field model also 
allows us to identify genuine MHD effects. 

Figure [7] shows the evolution of the central density (left 
panel) and the amplification of the magnetic energy (right 
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Fig. 8. Radial profiles of tlie angular velocity at the equator for model s20AlB5-D3M10 (left panel) and s20AlB5- 
D3M12 (right panel) at different times after bounce: t — tt, — 0, 8, 34, 45 and 59 ms. The time-independent rotation 
profile of the passive field model s20AlB5-D3M0 is shown by the black line in both panels. The yellow dashed line in the 
right panel indicates a change of sign of the angular velocity. 



panel) for all three models. The latter quantity is computed 
from 

En..^^\jd^x^Wh\ (43) 

Additionally the lower panel of Fig. [7] shows the L^c norm 
of (7b defined as the ratio of the total magnetic flux at 
the surface of each numerical cell to the average mag- 
netic flux on the surface. This dimensionless quantity mea- 
sures the quality of the numerical preservation of the diver- 
gence of the magnetic field along the evolution. The final 
value is consistent with the round-off error in the evolution, 
which can be computed as (double precision accuracy) x 
^(number of iterations) = lO^^^ x V5 x 10*^ = 2.3 x lO^^^ 
, if one considers a binomial distribution of errors. As the 
collapse proceeds both the density and the magnetic energy 
grow very similarly in all three models, because even in the 
highly magnetized progenitor model s20AlB5-D3M12 the 
strength of the magnetic field is insufhcient to affect the 
collapse dynamics. The rati o of magnetic energy to gravi- 
tational binding energy (see lCerda-Duran et al.|[2007l ) does 
not exceed a value of 10~^ (10"'^) during the collapse in 
model S20A1B5-D3M10 (s20AlB5-D3M12), which justifies 
the use of the passive field approximation in the weak mag- 
netic field limit. After core bounce, the low magnetized 
model s20AlB5-D3M10 continues to behave similarly to 
model s20AlB5-D3M0, since the magnetic field remains 
weak. The central density is slightly higher than in model 
s20AlB5-D3M0, but the magnetic field is far from satura- 
tion and is still growing linearly with time at the end of the 
simulation. 

On the other hand, the highly magnetized model 
s20AlB5-D3M12 clearly shows a saturation of the magnetic 
field energy shortly after core bounce. At this time the ra- 
tio of magnetic energy to gravitational binding energy is 
7%, a value that is never exceeded during the evolution. Its 
central density continues to grow beyond bounce, and the 
model eventually approaches an equilibrium configuration 
with a central density about 10% larger than in the passive 
field case. The behavior of the central density can be under- 
stood by examining the angular velocity profiles in Fig.O 
At the time of bounce, the angular velocity profile is very 
similar for all models, since the magnetic field is still unim- 
portant for the dynamics: the innermost 10 km of the core 



rotate rigidly, while further out follows a power law with 
an exponent ~ —1.2. This profile remains unaltered during 
the subsequent evolution of the passive field model. In the 
magnetized models, however, the central region spins down, 
and the central density rises, the effect being more promi- 
nent in the stronger magnetized model s20AlB5-D3M12. 
The right panel of Fig.[H] shows that the angular velocity 
begins to decrease for 10 km < r < 30 km shortly after 
bounce. In this region, the magnetic field is strongest since 
differen tial rotation winds up th e magnetic field more effi- 
ciently (jCerda-Duran et 311120071 ) . On a time scale of about 
50 ms, the angular velocity decreases by about a factor 10, 
and the innermost few kilometers of the core even acquire 
retrograde rotation. The reason for this effect is the increas- 
ing magnetic tension in the wound-up magnetic field lines. 
The characteristic time scale in which this magnetic ten- 
sion acts on the fluid is related to the Alfven crossing time 
scale of the innermost region ta ^ 50ms, which coincides 
with the time it takes for the retrograde rotation to appear. 
Th is effect was already o b served in N ewtonian simulations 
by iMiiller fc Hillebrandtl ()1979[ ). and lObergauhnger et al] 
()2006aD . For model s20AlB5-D3M10. the spin-down occurs 
more slowly, and saturates about 50 ms after bounce. 

To demonstrate the spin-down more clearly, we plot, 
in Fig. [21 the evolution of the central angular velocity for 
all three models. In the passive field approximation (black 
line), oscillates after bounce in accordance with the os- 
cillations of the core, and approaches a constant value at 
the end of the simulation. As the magnetic field increases 
in the progenitor, the spin-down of the core occurs more 
rapidly. This may be understood by means of the magneto- 
rotational instability (MRI hereafter). The MRI is a shear 
instability that can appear when both magnetic fields an d 
differential rotation are present (|Balbus fc HawlevI Il99l[ ). 
and it gives rise to transport of angular momentum. A 
necessary condition for the occurrence of the MRI is 
vjdrn^'^ < 0, where m — rsinO. In unstable regions, 
the MRI grows exponentially for all length scales larger 
than a critical length-scale Acrit ~ 27rcA/f^, where ca 
is the Alfven speed. The fastest-growing MRI mode de- 
velops on length-scales near Acrit on a typical time-scale 
of tmri = 4:7r[ujd^fi\~^ . Therefore, in order to numeri- 
cally capture the MRI, one has to resolve length-scales 
of about Acrit- Once the MRI grows, it develops chan- 
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Fig. 7. Evolution of the central density pc (upper panel), 
of the amplification of the magnetic energy i?mag/^'mag,o 
(middle panel) and the Loo norm of (Tb (lower panel) 
for model s20AlB5-D3M12 (solid), and s20AlB5-D3M10 
(dashed), respectively. The results obtained with the 
passive field appr o ximat ion (model s20AlB5-D3M0 of 
ICerda-Duran et al.l (|2007D ) are shown with a dotted line 
(almost overlapped to the dashed line in the middle panel). 



nel flows CHawl ev fc BalbuifTOOl. which are unstable to 
non-axisymmetric instabilities ( Goodman fc Xulll994D and 
eventually b ecome turbu lent in three-dimensional simula- 
tions (|Hawlev et al.lll996[ ). 

In our simulations, the region with r > 10 km is unsta- 
ble to the MRI due to its negative angular velocity gradi- 
ent. The growth time of the fastest-growing mode is in the 
range 1 ms to 10 ms for the region behind the shock wave, 
and about 1 s or even larger further outside. Since the time 
scale is independent of the initial magnetic field strength, 
these values are similar for both magnetizations (s20AlB5- 
D3M10 and s20AlB5-D3M12), and for the passive field 
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Fig. 9. Evolution of the central angular velocity Q.c for 
model S20A1B5-D3M12 (red), s20AlB5-D3M10 (blue), 
and model s20AlB5-D3M0 (black), respectively. The red 
dashed line (model s20AlB5-D3M12) indicates a change of 
sign of the angular velocity. 
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Fig. 10. Magnetic field topology at the end of the simula- 
tion, 51 ms after bounce, for model s20AlB5-D3M10 (upper 
panel) and s20AlB5-D3M12 (bottom panel), respectively. 
The ratio of magnetic to thermal pressure Pma.g/P is shown 
color-coded. Thin, white lines are poloidal magnetic field 
lines, and the thick, white line marks the neutrino-sphere. 
The axis labels are in units of km. 



case (s20AlB5-D3M0). However, the critical length scale 
depends on the strength of the magnetic field. 

For model s20AlB5-D3M12, the critical length scale at 
bounce is between Adit 1 km and 5 km inside the un- 
stable region (10 km < r < 30 km). This region is covered 
with 60 radial and 30 angular zones, which corresponds to 
a resolution {Ar,rA9) of 125 m x 500 m at r = 10 km, and 
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900 m X 1500 m at r = 30 km. This resolution is marginally 
sufficient to resolve the length scale of the fastest-growing 
mode of the MRI at bounce (5 — 10 radial zones, and 
2 — 3 angular zones). The strong redistribution of the angu- 
lar momentum observed for model s20AlB5-D3M12 might 
therefore be caused by the MRI. In turn, the saturation 
of the magnetic field is a direct consequence of this redis- 
tribution of the angular momentum. Without differential 
rotation, the poloidal magnetic field cannot be wound up 
into a toroidal magnetic field. The typical spin-down time 
scale Tspin_down Can be measured from Fig. [5] by fitting an 
exponential to the declining part of the curve. For model 
S20A1B5-D3M12, one obtains Tspin— down — 22.5 ms, which 
corresponds roughly to the time scale of the MRI. 

On the other hand, for model s20AlB5-D3M10 the crit- 
ical length-scale at bounce is about a factor of 100 shorter, 
i.e. between Acrit ^ 10 m and 50 m, and thus the fastest- 
growing mode of the MRI cannot be resolved with our grid 
resolution. Only modes with slower growth rates can be 
resolved on the grid. Accordingly, the spin-down for this 
model occurs on a longer time scale of Tspin-down = 62.9 ms. 
At about 50 ms after bounce, the innermost 10 km of the 
core develops a positive angular velocity gradient (see 
Fig-®, and hence becomes stable to the MRI. The cen- 
tral core is no longer able to lose angular momentum, and 
its spin down stops. We suspect that the appearance of this 
positive gradient is due to the poorly resolved MRI, which 
turns out to be more efficient in the inner region, where the 
resolution is higher, instead of where the shear is larger. 
The magnetic field continues to grow at similar rates until 
the end of the simulation due to the further winding-up of 
poloidal magnetic field lines, and because angular momen- 
tum transport is insufficient to affect the rotation profile 
outside the innermost 10 km significantly. 

Figure [TUl displays the magnetic field topology for mod- 
els S20A1B5-D3M10 and s20AlB5-D3M12 at the end of the 
simulation. The low magnetized model s20AlB5-D3M10 
(top panel) has a similar field structure as model s20AlB5- 
D3M0 of , Cerda-Duran et all 12007), since the passive field 
approximation holds very well for weakly magnetized pro- 
genitors (apart from its inability to capture the MRI). 
The prompt convection developing after bounce twists 
the magnetic field outside the neutrino-sphere, which is as- 
sumed to be located at pi, = 2 x lO^^gcm^-^, at about 
30 km. In model s20AlB5-D3M12, the magnetic field grows 
to values close to equipartition, and a distinctive, strongly 
magnetized outflow propagates along the axis behind the 
shock front. Between 10 km <^r < 30 km, where the MRI is 
predominantly growing, axisymmetric channel flows form, 
which are morp hologically similar to the flows found in the 
simulations of Hawlev fc Balbuj (119921 ). We analyze this 
issue in more detail in Fig. llli where the development of 
the channel flows is shown. Their length scale increases as 
the magnetic fleld becomes stronger, and since we assume 
axisymmetry they are stable, i.e. they do not cause any tur- 
bulence. 

Another important difference between models s20AlB5- 
D3M10 and s20AlB5-D3M12 is the location of the shock. 
At ^ 50 ms after core bounce, the shock is located 



^ This transient is produced by an unstable entropy gradient, 
which is probably an artifact of our poor neutrino treatment . 
The interested reader is addressed to lCerda-Duran et al.l (|2007f ) 
for a detailed discussion of this issue. 



about 50 km further out in the strongly magnetized model 
s20AlB5-D3M12 than in the weakly magnetized model 
s20AlB5-D3M10. This is most likely a consequence of the 
transport of angular momentum by the MRI which pushes 
the shock front to a larger radial distance, although our cur- 
rent grid resolution is probably too poor in the shock region 
to conflrm this interpretation conclusively. Understanding 
this effect and, in particular, its implications for the explo- 
sion mechanism, requires a separate study, which will be 
published elsewhere. 

Since no other simulations yet have been published that 
are capable of treating a similar combination of general rel- 
ativity and microphysics, it is impossible to compare di- 
rectly with other work. Nevertheless, we find qualitative 
agreement with related s imulations of ma g neto-rotational 
core collapse dObc rgaulin eer et al.l l2006bl [al: IShibata et al.l 
2001 II 

urrows et al.l l2007). In particular, our simulations 
share the following aspects with these investigations: (i) re- 
distribution and transport of angular momentum radially 
outwards due to the MRI, resulting in the spin down of 
the central region of the core; (ii) increase of the central 
density after core bounce due to angular momentum losses; 
and (iii) appearance of a weakly relativistic but highly mag- 
netized outflow along the axis. This agreement strengthens 
our confidence in the suitability of our new numerical code 
for the systematic investigation of magneto-rotational core 
collapse, which we shall report elsewhere. 

5. Conclusions 

We have presented a new numerical code that solves the 
GRMHD equations coupled to the Einstein equations for 
the evolution of a dynamic spacetime. Hence, it extends the 
small list of available codes that are capable of modeling 
these challenging physics. The main objective of the new 
code is the study of astrophysical scenarios in which both 
strong magnetic fields and strong gravitational fields are 
present, such as the magneto-rotational collapse of stellar 
cores, the coUapsar model of GRBs, and the evolution of 
neutron stars. 

Our new numerical code is based on high-resolution 
shock-capturing schemes to solve the flux-conservative hy- 
perbolic GRMHD equations, and the constraint-transport 
method to ensure the solenoidal condition of the magnetic 
field. The Einstein equations are formulated in the CFC ap- 
proximation, and the resulting elliptic equations are solved 
using a linear Poisson solver. The motivation to use CFC 
is based on the astrophysical applications envisaged for the 
code, which do not deviate significantly from spherical sym- 
metry. Furthermore, the code incorporates several equa- 
tions of state, ranging from simple analytical expressions 
to tabulated microphysical equations of state. 

We have presented a number of stringent tests of our 
new GRMHD numerical code, which are the main focus 
of this paper. The test calculations demonstrate the abil- 
ity of the code to handle properly all aspects appearing in 
the astrophysical scenarios the code is intended for, namely 
relativistic shocks, strongly magnetized fluids, and equilib- 
rium conflgurations of magnetized neutron stars. One of the 
tests the code has passed successfully is in fact an applica- 
tion, namely the simulation of general relativistic magneto- 
rotational core collapse using a realistic stellar progenitor 
model and a microphysical equation of state. We have com- 
pared the results obtained by our new code with those of 
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Fig. 11. Details of the magnetic field structure of the core at three different times after bounce: t — ti, — 9 ms (left), 
11.5ms (middle), and 14ms (right), respectively. The ratio of magnetic pressure to thermal pressure Pmag/P is shown 
color-coded. Thin, white lines are poloidal magnetic field lines, while the thick, white line marks the neutrino-sphere. 
The axis labels are in units of km. 



a previous study based on the passive magnetic-field ap- 
proximation, and find good agreement for initially weakly 
magnetized progenitors. 

Finally, we mention that the new code is also capable of 
handling the gravitational collapse leading to the formation 
of a black hole. Results for this specific application will be 
presented elsewhere. Further extensions of the code that 
we foresee in the near future include the incorporation of 
a simplified scheme for neutrino transport (to explore the 
post-bounce evolution of collapsing magnetized cores more 
reliably) along with the implementation of resistive MHD. 
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